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ABSTRACT 

We study the two dimensional Kolmogorov flows in the limit as the forcing frequency 
goes to infinity. Direct numerical simulation indicates that in this limit the low frequency 
energy spectrum evolves to a universal k~ 4 decay law. We derive effective equations gov- 
erning the behavior of the large scale flow quantities. We then present numerical evidence 
that with smooth initial data, the solution to the effective equation develops a k~ 4 type 
singularity at a finite time. This gives a convenient explanation for the k~ 4 decay law 
exhibited by the original Kolmogorov flows. 
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§ 1. Introduction 


Toward the end of the 50’s, Kolmogorov proposed to study the two-dimensional incompress- 
ible flow with a high frequency forcing term as a model to help understand the transition to 
turbulence and the inverse cascade process in a turbulent flow. These so-called Kolmogorov 
flows are described by the following equation: 


( 1 . 1 ) 


ut + u • Vu + Vp = v/\u + vk 2 j 

< 


sm(kfX2 ) 
0 


) 


V • u = 0. 


Here u is the dynamic viscosity, k/ is the forcing frequency. Since then, this problem has 
attracted a lot of attention. Considerable progress has been made on both the experimental 
and theoretical side. We refer to [10] for a discussion on the laboratory realization of 
the Kolmogorov flows. The theoretical side of this problem has mostly been concentrated 
on the hydrodynamic stability viewpoint. Linear stability of the basic Kolmogorov flow 
(a stationary shear flow proportional to the forcing) was studied in [6]. The amplitude 
equations for the weakly nonlinear theory were derived in [9, 15]. It was found that this 
problem belongs to a wider class of problems which exhibit large scale instability. Later stage 
development of the Kolmogorov flow was extensively studied more recently in [11, 12] using 
a dynamical system viewpoint. Several scenarios of the transition to chaos and turbulence 
were identified by resorting to careful numerical simulations. 

In this paper, we take a different approach and study Kolmogorov’s problem in the 
following asymptotic limit: 

v = y~ = £ — »• 0 . 

k f 

Our main interest is to understand the inverse cascade process in the Kolmogorov flow. Our 
extensive numerical results indicate that, even if we start with a smooth initial data with 
an exponentially decaying energy spectrum, the high frequency forcing term pumps energy 
into lower frequencies so that, beginning at a finite time, the energy spectrum exhibits a 
universal k~ 4 algebraic decay law at the low frequencies. 
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To explain this phenomenon, we derive the effective equations governing the evolution of 
the large scale flow quantities. We present strong numerical evidence that the solution of the 
effective equation with smooth initial data form a singularity at a finite time. Furthermore, 
the energy spectrum at the time of singularity formation is precisely k~ 4 . 

We remark that in this paper the term “cascade” is used loosely to describe situations 
where either energy or enstrophy is shifted in the frequency space. We do not have evidence 
that energy or enstrophy is transferred continuously in the frequency space with nonzero 
flux at each wavenumber since we did not measure these quantities. 

This paper is organized as follows. In the next section we present results of direct 
numerical simulation of the Kolmogorov flow in the limit as the forcing frequency goes 
to infinity. The emphasis will be placed on the evolution of the energy spectrum. In 
Section 3, we derive the effecetive equations for the averaged quantities using multiple scale 
asymptotics. In Section 4, we present our numerical evidence for the singularity formation in 
the solutions of the effective equations. We conclude this paper in Section 5 by commenting 
on the relevance of our findings for more general flows, including three-dimensional flows. 

§ 2. The Inverse Cascade Phenomenon in the Kolmogorov Flows 


As we mentioned in the introduction, we are interested in studying the inverse cascade 
phenomena associated with the Kolmogorov flow in the following asymptotic limit: 


(2.1) v = — = e. 

k f 

As can be easily seen from the derivations presented in the next section (see Remark 3.2), 
this is the only distinguished limit for which both the viscous force and the external high 
frequency force play important roles in the effective equations governing the evolution of 
large scale structures. 

Under the scaling (2.1), equation (1.1) takes the following form: 


( 2 . 2 ) 


+ u e • Vu £ + Vp e = eA u £ + 


1 

£ 


( sinf- ^ 

0 / 


V • u s = 0 
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To study the influence of the small-scale forcing on the evolution of the large scales, in 
particular, to study the evolution of the energy spectrum associated with the Kolmogorov 
flow, we have numerically integrated (2.2) using a number of initial velocity fields and 
studied the limit as e —*■ 0. We discretize (2.2) using a Fourier-collocation method in space, 
and a third order Runge-Kutta method in time. More details of the numerical method will 
be given in Section 4. Our calculation differs from the previous calculations of Sulem et.al. 
[16] and Platt et.al [11] in that we studied as far as we can the limit as e — ► 0. In particular we 
computed the solutions of (2.2) with forcing wavenumbers kj — 10, 25, 50, 100 (corresponds 
to e = 0.1,0.04,0.02,0.01). In these cases the solutions exhibit a clear separation of scales 
for some time. 

Here we report two examples of our computational results. We have experimented with 
other sets of initial data. In all cases we observe the same kind of phenomena as we report 
below. 

In the first example we take the initial vorticity to be 


(2.3) 


du 0 dvo . 2 x l • 2 x 2 ? 

Wo = -5 Q — = sin 2 sm 2 - 7r 2 

dx 2 dxi 2 2 


We computed the energy spectrum using the following formula: 


(2.4) E(k) = £ (K*i,* 2 )| 2 + K*i,*2)| 2 ) 

where u(fci,& 2 ) is the (&!, fc 2 )-th Fourier mode of u. Figure 1 is the log-log plot the energy 
spectrum at times tj = 0.1, 0.2, * ' * , 2.5. Here e = 0.04, and the computation is done on a 
400 2 grid. It is clear that the energy spectrum consists of two parts: The high frequency 
part is dominated by the forcing frequency and its higher harmonics. The low frequency 
part corresponds to the large scale structure of the flow. Although it is not shown in Figure 
1, the energy spectrum at t = 0 decays much faster than those plotted there. At early times 
the low frequency spectrum is very well separated from the high frequency part. At about 
t = 1, the high end of the low frequency spectrum begins to merge with the low end of the 
high frequency spectrum. At the same time, the low end of the low frequency spectrum 
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forms an envelope. This envelope is compared to a straight line with slope equal to —4. 
They are close. 

In the second example we choose a nonsmooth initial vorticity: 

(2.5) u(x\,x-i) = H(X\)H(X2) - 7T 4 , 

where H(x) is the hat function on [0,2x] with height 7r. The other parameters are the same 
as in the first example. The energy spectrum for t = 0.1, 0.2, • • • , 1 are plotted in Figure 2. 
We observe basically the same phenomena as in Figure 1 except that here the initial data is 
much less smooth, therefore the low frequency spectrum settles down to an envelope much 
faster. The envelope is also compared to a straight line of slope -4 and is seen to be very 
close to that. 

Although there exists a large literature on two-dimensional turbulence theory (for a 
review, see [2]), the discussions are much less sound than Kolmogorov’s theory for three- 
dimensional turbulence. We therefore seek our own explanation of this seemingly universal 
behavior of relaxation to a k~ 4 law for the low frequency spectrum. In the next section, we 
will derive the effective equations governing the evolution of the large scale structures and 
therefore the low frequency part of the spectrum. In Section 4 we will study the effective 
equations numerically and present evidence that the energy spectrum for the solutions of 
the effective equations evolve to a k~ 4 decay law, even with a smooth initial data. The time 
that this happens is close to the time when the envelope forms at the low frequency part of 
the spectrum of the Kolmogorov flows. Beyond this time the original Kolmogorov flows no 
longer exhibit separation of scales and the effective equations cease to be valid. 

§ 3. Effective Equations for the Large Scales 

In this section, we derive effective equations governing the evolution of large scale flow 
quantities. Enormous efforts have been devoted to the the derivation of similar equations, 
the so-called “Reynolds equation”, for general turbulent flows. The fact that turbulent flow 
has a continuum range of scales presents a serious difficulty to this objective. Our problem 
is drastically simplified since there is a separation of the scales which are active in the flow. 
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As we pointed out in the last section, eventually separation of scales no longer holds and 
the flow exhibits a transition to a continuum of scales. This signals the breakdown of the 
effective equation, an issue which will be looked at more carefully in the next section. 

We begin with the following proposition: 


u £ (x,t) 

v s (x, t) 


u(x, t) 
v(x ,t) 


?(x,t, 


7 w(x,t,^) j + £ | 


V t (x,t,f) 


p e (x, t ) = p(x, t ) + II(a;, t, — ) + £ Pl (x, <,—) + ••• 

£ £ 

where w(x,t,y 2 ), U(x,t,y 2 ), ui(x,t,y 2 ), p\(x,t,y 2 ) etc. are periodic functions of y 2 , and 
< w(®> V) >= o, < z(x,t,-)>= 0, < n(®,<,-) >= 0 for every *, t 6 R 2 xR + . Substituting 
(3.1) into (2.2), and collecting equal orders of e, we get a hierarchy of equations. The 0(e~ l ) 


equation are: 


(3.2) V, 


(u + w) 2 (u + in)(u + z) 
( u + w)(v + z) (v + w ) 2 


+ Vj/II = Aj, 


w / sin 2/2 

+ 


Vy(w,z) T = 0. 

The solutions of (3.2) can be obtained explicitly: 


w = ^ 2~T { [ sin V 2 ~v cos y 2 ], z = 0, II = 0. 


We next come to the 0(1) equations: 


u(x,t) I w ' 

+ +V X 

v(x,t) ) z / ^ 


(u + w ) 2 (m + «>)(« + z) 

(• u + w)(u + z) (h + in) 2 


( 2(u + w)u i (u + w)v i + (v + z)u\ \ 

v„- + v K p + v^n + v y px 

\ (u + w)vi + (v + z)u\ 2(v + z)v i ) 


W Ml 

= 2(V. • V„) ■ ■ ■ + A v 


Vj. • (u + w, V + z) T + Vy • (Ml, V\) T = 0. 
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Averaging both sides of (3.4) with respect to y, we get (written in conservative form): 


(3.5) 



/ <-2 .. 

' u i uv 


UV V 


\ 


) 


+ v a 


< w 2 > <wz> 
<wz> < z 1 > 


\ 


+ V x p = 0 


Substituting in (3.5) the result of (3.3), we get 

f 1 

Ut + UU X1 + vu X2 + ( ^T+l)) 371 + P* 1 ~ 0 
( 3 - 6 ) ' v t + uv Xl + VV X2 + p X2 = 0 

U x J + v X2 = 0 

This is the effective equation governing the evolution of the large scale quantities for the 
Kolmogorov flow. For the self-consistency of this asymptotic argument, we have to produce 
a candidate («i,pi) which solves the remainder of (3.4). To this end we write (3.4) as 


/ 


(3.7) V t 

(3.8) 
where 

-/ = 


2(u + w)u\ ( it + w)v i + (v + z)ui 

(it + w)v\ + (D + z)u\ 2(f) + z)v i 

Vj ) • («i, vi) T = -V x • (u + w,v + z) T = g 


4 " ^ yPl 



-a] 

-h) 


( u(x,t ) ‘ 



\ KM) j t 

+ + VJI — 2(V* • Vy) 


( u + w ) 2 
^ (ii 4- w)(v + z) 


(u + w)( V + z) 
(f) + w ) 2 




In (3.7) and (3.8), a; and t are viewed as parameters. Since / and g only depend on j/ 2 , 
we seek a solution of (3.7) and (3.8) which depends only on y 2 - The constructions above 
for w, z, II, u, v and p guarantee that the averages of / and g in t /2 over the period [0,27 t] 
are zero. Therefore we solve v\ from (3.8). The first equation of (3.7) will then give the 
solution U\ whereas the second equation will give the solution pi . It is easy to see that if 
we define 


(3.9) 



^ u(x,t) ^ 

KM) } 


+ 


KM,^) ) 


+ £ 
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p s (x, t) = p(x, t) + n(®, t, + epi(x, t, ^), 

then the equations (2.2) are satisfied by u e , v s ,p e with an error of order e. 

Remark 1: The difference between equation (3.6) and the usual 2-D Euler’s equation 

is the appearance of the off-diagonal term in (3.6). This means that the linearized equation 
is weakly well-posed: high frequency perturbations will grow with a growth rate which is 
a polynomial in the wavenumber. This fact motivates us to look for finite time singularity 
formation in the solutions of (3.6) with smooth initial data. 

Remark 2: If we change the forcing term in (1.1) to a different function, the form 

of the effective equation will not change. The off-diagonal term will usually change to a 
different function. In particular if the viscosity is scaled to vanish faster than 0(e), then 
the added flux in the effective equation is singular in v. 

Remark 3: We notice that (3.6) conserves energy: 

\ d f f 1 

2dt J + dXldX 2 = ~J ^ 2(fi 2 + ' i )^ dxidx2 

= / *- 2tvhv ) iXliX2 = - J WthT ) iXliX2 = °- 

§ 4. Singularity Formation in the Solutions of the Effective Equations 

In this section, we present numerical evidence which suggests that the solutions of (3.6) 
with analytic initial data develop singularities at finite time. Furthermore, at the time when 
singularity forms, the energy spectrum associated with the solutions decays like k~ 4 . 

The numerical results presented below were computed using a Fourier-collocation method. 
We checked these results using both a second order centered differencing scheme and a fourth 
order ENO type scheme [14]. We remark that searching singularities numerically is a very 
subtle undertaking. It is important to check the numerical results using different meth- 
ods since each method is likely to have its own numerical artifact when it comes to the 
computation of singularities. 

Strictly speaking, demonstrating finite time singularity formation is an analytical prob- 
lem. The ultimate answer has to come analytically, either by giving specific examples or 
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by establishing a lower bound for the life span of the smooth solutions. Computers are 
limited to finite arithmetic and finite capacities. A numerical approach can at best give 
partial evidence, not the complete answer. Needless to say that we have tried the analytical 
approach and the problem seems to be very hard. For the moment, presenting numerical 
evidence is the best we can do to shed some light on our problem. 

§ 4.1. Description of the Numerical Method 

For the spatial discretization, we used the Fourier-collocation method [1], Roughly 
speaking, the differentiation operator is approximated in the Fourier space, whereas the 
nonlinear operations such as multiplication are done in the physical space. We use the 
intrinsic Cray FFT routines which considerably enhanced the performance of the code. 

We observed that even for the computation of smooth solutions of the standard two 
dimensional incompressible Euler’s equation, a certain amount of filtering is needed for the 
stability of the numerical method, and it is crucial to add the filters correctly. A robust way 
of adding the filters is to replace the Fourier multiplier ikj for the differentiation operator. 
■for by ikj<p(\kj\), where 

(4.1) <p(k) = e~ a( $ r) " V , for \k\<N 

Here N is the numerical cutoff for the Fourier modes, mf is the order of the filter, a is 
chosen so that <p(N) = e~ a = machine accuracy. The machine accuracy on Cray-YMP with 
single precision is roughly 10 -14 . The accuracy of such an approximation scheme depends 
on the parameter m,f. For smooth functions f(x), we have 

(4.2) \\f'{x)-D N f{x)\\ = 0(N- m s) 

With -Djv/ = (ikip{\k\))F f is the numerical approximation of f'{x). F denotes the 

Fourier transform operator. Unless otherwise stated, the results presented in § 4.2 were 
computed with m/ = 14. 

For the temporal discretization, we used the third order Runge-Kutta method designed 
in [14], We used the third order version since it only requires three auxiliary arrays, whereas 
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the fourth order version requires five auxiliary arrays. We take initial data that is periodic 
with period D where D = [0,2tt] x [0,27t]. 

§ 4.2. Numerical Results 

We present our numerical results for the case when the initial vorticity is the same as 
in the first example of Section 2: u 0 = sin 2 ^ sin 2 | - 7 r 2 . 

As a first hint for the singularity formation, we plotted the time history of the quantity 
H\(t) — J D \Vu\ 2 dxdy. Figure 3 suggests that as N —*■ +oo, Hi(t) grows without bound, 
for t > 0.9. 

This information alone cannot be used as a solid evidence for singularity formation, 
since it only says that the quantity Hi(t) is poorly resolved for t > 0.9. Next we plotted 
in Figure 4 the time evolution of the energy spectrum (on a log-log scale) associated with 
the numerical solutions obtained with three different resolutions: 128 2 , 256 2 and 400 2 . We 
observe that at t = 1, a bump develops near the cutoff frequency on the energy spectrum. 
This indicates that at this time, the numerical solutions are not very well resolved. We also 
observe that the energy spectrum at t = 0.9 is well represented by the numerical solutions 
obtained on the 128 2 and 256 2 grids, and it fits very well with a straight line with slope 
—4. This is a stronger indication that at this time, the energy spectrum of the solution no 
longer decays exponentially, but decays algebraically as A: -4 . 

To determine more precisely when algebraic decay sets in, we use the techniques pro- 
posed by Sulem, Sulem and Frisch [17]. Let f(z) be a function which is analytic in the strip 
|Irn2| < a, and has a branch-cut singularity z 13 on |Imz| = a, then its Fourier transform 
f(k) decays like 

(4.3) f(k) ~ Ck~ 0 e~ ak , |*| -► +oo. 

Imagine the following scenario for the singularity formation in the solutions of (3.6). 
Instead of solving (3.6) in the real space R 2 , we solve (3.6) in the complex space ( 21 , 22 ) £ C 2 . 
Initially the solution is an entire function. At later time the solution develops branch point 
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of order /3(t ) at z(t ) with |Im^r(t)| = a(t) > 0, but it is analytic in the strip |Im«| < a(f). 
If this branch point travels to the real axis at a finite time t*, then the solution of (3.6), 
solved in the real space, develops a singularity at t = t*. 

This scenario has been proposed for a number of problem in fluid mechanics, notably 
the formation of singularity in the evolution of a vortex sheet, (see [8, 3, 13]). In this case, 
there are strong numerical as well as analytical evidence which confirms this picture. 

To check whether this provides a plausible picture for our problem, we checked numer- 
ically the validity of proposition 4.3 for the energy spectrum. To do that we choose three 
consecutive wave numbers ( k - 1, k,k + 1), and solve 

(4.4) E(k,t) = Ck-0e~ ai t for k = k - 1, k, k + 1, 

with double precision arithmetic (on Cray-YMP). At each fixed time, we get three functions 
0(k,t), a(k,t). These functions depend on k. In order that (4.3) be a good estimate 
for the energy spectrum of the solutions of (3.6), asymptotically as k -+ oo, the functions 
C(k,t), (3(k,t ) and a(k,t ) should not depend on k in this limit. In Figure 5 we plotted the 
functions a(k,t ) from t = 0.4 to t = 1 with constant time increments. The computation is 
done on a 256 2 grid. It is clear from Figure 5 that for fixed t, a(k,t) is almost a constant for 
20 < k < 60. For k < 20, a(k, t) exhibits fluctuations since k is not yet in the asymptotic 
regime. For 80 < k < 128, the energy density is at the order of 10 -32 which is below machine 
accuracy even if double precision is used. These wavenumbers are needed only to enforce 
smooth decay for the Fourier coefficients down to the machine accuracy. Figure 6 displays 
similar results for a computation using a 400 2 grid. The fluctuations at high wavenumber 
is a manifestation of trying to fit machine zeros with proposition 4.3. The functions C(k,t) 
and /?(&, f ) exhibit the same feature. 

Now we can average the values of a over the intermediate wave numbers and safely 
regard the averaged value (as a function of t ) as the width of the analyticity strip of the 
solution at that time. This function will be denoted by a(t). The corresponding function 
for /? will be denoted by /3(f). Displayed in Figure 7 is the function a(t) for two different 
numerical resolutions: 256 2 and 400 2 . Results obtained on both grids indicate that the 
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width of the analyticity strip for the solution vanishes at t* ~ 0.9, suggesting a singularity 
formation at this time. The corresponding results for (3 are displayed in Figure 8. As t -+ t*, 
(3(t) converges to 4. 

The other numerical parameter in our method is the order of the filter. Figures 9 and 10 
display the results for a, /? as we vary the order of the filter. We see that the approximated 
values of a, f3 are not sensitive to the change of the filters. However when the filter is too 
weak (my = 15), the numerical solution deteriorates drastically at a time before /■*, since 
it can no longer handle the large gradients developed in the vorticity field. The computed 
value fluctuates much more, therefore the result is not shown in Figure 10. This is not 
unexpected since in general /3 is much more sensitive than a. 

To get an idea about the nature of the singularity, we display in Figures 11 and 12 the 
time evolution of u and u at y — ^ with t = 0.3, 0.6, 0.9. At t — 0.9, a cusp seems to have 
formed on the profile of u>. Contour plots of u show that these cuspy structures occur on 
isolated lines. 

The behavior of the solutions reported above seems to be generic. It occurs in other 
calculations we did using different initial data. 

§ 5. Concluding Remarks 

As we mentioned in the introduction, the Kolmogorov flow belongs to a large class of 
problems which exhibit large scale instability. A typical example in three dimension is 
the ABC flow which is known to be unstable to large scale perturbations [7]. In principle 
the ideas and methods presented above should apply to these problems also, although to 
actually carry out the program is difficult mainly because the computational cost increases 
drastically for three-dimensional problems. However, it is clear that such instabilities will 
trigger the transport of energy from small scales to large scales, thereby effecting the mean 
flow quantities. 

One can also attempt to study the interaction of different scales by studying the solutions 
with two-scale initial data. Formally one can use multiple scale asymptotics to derive 
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effective equations for the large scales. This is done in [5, 4]. However, a closer look at 
these effective equations reveals that the linearized equations are exponentially ill-posed: 
high frequency perturbations grow at a rate that is exponential in the wavenumber. The 
reason for this is very simple. The incompressibility condition forces the microstructure 
to be shear flows. In the limit of infinite scale separation, the microstructure appears like 
stacks of vortex sheets when viewed at macro-scales. The linear ill-posedness is simply a 
manifestation of the Kelvin-Helmholtz instability. 

In such situations multiple scale asymptotics does not provide useful information since 
it is based on the hypothesis of scale separation. For the Kolmogorov flow, we are saved 
since there is a fixed time interval in which the scales are separated and we can get useful 
information from two-scale expansions. As one might expect, beyond the critical time t *, the 
two scales in the Kolmogorov flow merge and the effective equation cease to be valid. This is 
checked numerically with the following procedure. The effective equation (3.6) is integrated 
numerically using a finite difference method with a built-in nonlinear numerical viscosity 
to avoid oscillations coming from large vorticity gradient [14]. The energy spectrum of the 
numerical solution is found to be decaying increasingly slower. This differs sharply from the 
energy spectrum of the original Kolmogorov flow since the latter approaches an asymptote 
after the critical time. 
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Figure 1. Log-log plot of the time evolution of the energy spectrum with initial data (2.3) 
at t = 0.1, 0.2, •• • , 2.5. e = 0.04. t increases upward. The energy spectrum is compared to 
a line of slope -4. 


Figure 2. Log-log plot of the time evolution of the energy spectrum with initial data (2.5) 
at t = 0.1, 0.2, * • * , 1. £ = 0.04. t increases upward. The energy spectrum is compared to a 
line of slope —4. 
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Figure 5. Functions a(k, t ) from t — 0.4 to t = 1 computed from (4.4). t increases downward. 
Computation is done on a 256 2 grid. 



Figure 6. Functions a(k, t ) from t = 0.4 to t = 1 computed from (4.4). t increases downward. 
Computation is done on a 400 2 grid. 
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Figure 8. fi(t) computed on different grids. The line with circles is computed on 400 2 grid. 
The other line is computed on 256 2 grid. 
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Figure 9. The analyticity strip width a(t) computed on a 256 2 grid with different orders of j 
filter. The line with circles is computed with to/ = 5. The line with “x” is computed with 
to/ = 15. The other line is computed with to/ = 10. 



Figure 10. computed on a 256 2 grid with different orders of filter. The line with circles 
is computed with to/ = 5. The other line is computed with to/ = 10. The result for 
to/ = 15 fluctuates and is not shown here. 
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